cap program drop arsimsens
//version 1.2--3/17/15.


program define arsimsens,rclass
version 9

syntax varname [if] [in], Treat(varname numeric) Str(varname numeric) [Gam(real 1.0) Del(real 1.0) Verbose Hlestimate] 



if `c(changed)' == 1{
di "{error:Save data in memory before proceeding (note: clear option invalid for arsimsens)}"
exit 4
}

if `gam'<1{
di "{error:Enter a value of Gamma that is at least 1}"
exit 198
}

if `del'<1{
di "{error:Enter a value of Delta that is at least 1}"
exit 198
}

local mainfile= "`c(filename)'"

//begin quietly
quietly{


marksample touse
markout `touse' `treat' `str'

//for clarity in code:
local resp="`varlist'"

//convert Gamma to gamma, Delta to delta:
local udel=ln(`del')
local ugam=ln(`gam')

//for use of hlestimate only
//--------------------------------
if "`hlestimate'"=="hlestimate"{
summ `resp' if `touse'==1 & `treat'==1
local maxtreat=`r(max)'
local mintreat=`r(min)'

summ `resp' if `touse'==1 & `treat'==0
local maxctrl=`r(max)'
local minctrl=`r(min)'

local tau_l0=`mintreat'-`maxctrl'
local tau_r0=`maxtreat'-`minctrl'
}
//----------------------------------

//for aligned-rank test:
//stratum avg resp:
egen double __stavresp=mean(`resp') if `touse'==1, by(`str')
//"aligned response":
gen double __alresp=`resp'- __stavresp if `touse'==1
//aligned rank:
egen double __respvec=rank(__alresp) if `touse'==1
//gkr2000 adjustment:
quietly summ __respvec if `touse'==1
replace __respvec=__respvec/(.5*`r(max)') if `touse'==1

//observed test stat:
egen double __tv=total(__respvec) if `treat'==1 & `touse'==1
egen double __tv2=mode(__tv)
local tsobs=__tv2
drop __tv __tv2

//generate Z and r_C
gen double __rvec=__respvec
gen double __zvec=`treat' if `touse'==1


//generate the following vectors--to be used to calculate denominators of
//mu_ik and sigma^2_ik
//exp(udel*r_Ci):
gen double __edwvec=exp(__respvec*`udel') //touse implicit
//exp(ugam*Z_i):
gen double __egbvec =exp(`treat'*`ugam') if `touse'==1


 

sort __rvec


tempfile tfile
save "`tfile'",replace

quietly summ `str' if `touse'==1
local strmax=r(max)
local strmin=r(min)


//initialize running sum (over i) of mu_imax and sigma^2_imax:
tempname summu_imax sumsigma2_imax
scalar `summu_imax'=0
scalar `sumsigma2_imax'=0

//start loop 1: over i
forvalues i=`strmin'/`strmax'{ 
use "`tfile'",clear
marksample touse
keep  if `str'==`i' & `touse'==1

if `c(N)'==0{
continue
}

if "`verbose'"=="verbose" { 
noisily di "STRATUM: `i'" 
}

//`c(N)' works because only `touse' observations were kept for stratum i
//and c(N) is number of observations currently  in ds (even if not saved)
local n_i=`c(N)'



//this is to make factorials more readable in code:
//e.g., " 4! " is written as: " `fac'4`_' "
local fac="round(exp(lnfactorial("
local _=")),1)"


//calculate denominators for each mu_i`k', `k'=0...n_i:
//den 1:
//elsympolycalc calculates elementary symmetric polynomials e_k(Y_i)
//for vector Y_i --see below for program 
elsympolycalc __egbvec if `touse'==1
forvalues k=0/`n_i'{
local NmK=`n_i'-`k'
tempname den1_`k'
scalar `den1_`k''=(`fac'`k'`_')*(`fac'`NmK'`_')*(`r(e_`k')')
//thus den1_`k' is `k'!(n-`k')!*e_`k'(exp(ugam*Z_i))
}

//den2: 
elsympolycalc __edwvec if `touse'==1
forvalues k=0/`n_i'{
local NmK=`n_i'-`k'
tempname den2_`k'
scalar `den2_`k''=(`fac'`k'`_')*(`fac'`NmK'`_')*(`r(e_`k')')
//thus den2_`k' is `k'!(n-`k')!*e_`k'(exp(udel*r_Ci))
}


**************start mu_ik specific calcs***************************
//initialize running sum of inner sum Y_i
forvalues k=0/`n_i'{
tempname mu_i`k'
scalar `mu_i`k''=0 
}

summ __zvec
if `r(sum)'==1{ //i.e., if in stratum there is one treated obs
 local unique=1
 }
 else{ //i.e., if there is more than one treated obs --> one control obs
 local unique=0
 }
 
 

 gen double __zperm=abs(1-`unique') //i.e., a vector of 0s if one treated obs
 //and a vector of 1s if one control obs.
 
 //now cycle through vecs in orb(Z) 
forvalues x = 1/`n_i'{ //start permutation of Z
 local Xm1=`x'-1
 replace __zperm= `unique' in `x'
 capture replace __zperm=abs(1-`unique') in `Xm1'
 
 
 //calcuate a_tau--i.e., running sum of rho(Z)^TR
 tempname a_tau
 scalar `a_tau'=0
 forvalues xx=1/`n_i'{
  scalar `a_tau'=__zperm[`xx']*__rvec[`xx']+`a_tau'
 } // a_tau =rho(Z)^T*r when loop exits.
 
 //generate running sum  of inner sum:
 gen double __Y=exp(`ugam'*__zperm+`udel'*__rvec)
 elsympolycalc __Y
 forvalues k=0/`n_i'{
  local NmK=`n_i'-`k'
  scalar `mu_i`k''=`a_tau'*(`fac'`k'`_')*(`fac'`NmK'`_')*(`r(e_`k')')+`mu_i`k'' 
  //mu_i`k' will be mu_ik once next loop (over Orb(Z)) exits
  //i.e., once divided by denominators and multiplied by (n-1)!
 }
 drop __Y
} //end permutation of Z.
drop __zperm 
forvalues k=0/`n_i'{
local Nm1=`n_i'-1
scalar `mu_i`k'' = (`fac'`Nm1'`_')*(`mu_i`k'')/(`den1_`k''*`den2_`k'')
}

**************start sigma2_ik specific calcs***********************
//this parallels mu_ik calculation except term outside inner sum
//  is (rho(Z)Tr-mu_i`k')^2
//a little redundancy is present in these calculations in interest of clearer code.

forvalues k=0/`n_i'{
tempname sigma2_i`k'
scalar `sigma2_i`k''=0 
}

summ __zvec
if `r(sum)'==1{ //i.e., if in stratum there is one treated obs
 local unique=1
 }
 else{ //i.e., if there is more than one treated obs --> one control obs
 local unique=0
 }
 

gen double __zperm=abs(1-`unique') //i.e., a vector of 0s if one treated obs
 //and a vector of 1s if one control obs.


//now cycle through vecs in orb(Z) 
forvalues x = 1/`n_i'{ //start permutation of Z
 local Xm1=`x'-1
 replace __zperm= `unique' in `x'
 capture replace __zperm=abs(1-`unique') in `Xm1'

 
 //calcuate a_tau=rho(Z)^Tr
 tempname a_tau
 scalar `a_tau'=0
 forvalues xx=1/`n_i'{
  scalar `a_tau'=__zperm[`xx']*__rvec[`xx']+`a_tau'
 } // a_tau is pi(Z)^T*r when loop exits.
 //generate running sum  of inner sum:
 gen double __Y=exp(`ugam'*__zperm+`udel'*__rvec)
 elsympolycalc __Y
 forvalues k=0/`n_i'{
  local NmK=`n_i'-`k'
  scalar `sigma2_i`k''=(`a_tau'-`mu_i`k'')^2*(`fac'`k'`_')*(`fac'`NmK'`_')*(`r(e_`k')')+`sigma2_i`k'' 
  //sigma2_i`k' will be sigma2_ik once next loop (over Orb(Z)) exits. 
  //i.e. once divided by denominators and multiplied by (n-1)!
 }
 drop __Y
} //end permutation of Z.
drop __zperm
forvalues k=0/`n_i'{
scalar `sigma2_i`k'' = (`fac'`Nm1'`_')*(`sigma2_i`k'')/(`den1_`k''*`den2_`k'')
}



//***********calculate mu_imax**************
//*********..and sigma2_imax**************
tempname mu_imax sigma2_imax
scalar `mu_imax'=0 
scalar `sigma2_imax'=0
forvalues k=0/`n_i'{
if `mu_i`k''>`mu_imax'{
scalar `mu_imax'=`mu_i`k''
scalar `sigma2_imax'=`sigma2_i`k''
 }

if `mu_i`k''==`mu_imax' & `sigma2_i`k''>`sigma2_imax'{
scalar `mu_imax'=`mu_i`k''
scalar `sigma2_imax'=`sigma2_i`k''
 }
}

//running sum (over i) of mu_imax and sigma^2_imax: 
//(once exits loop over i, gives quantities used in (6) of Small et al)
scalar `summu_imax'=`summu_imax'+`mu_imax'
scalar `sumsigma2_imax'=`sumsigma2_imax'+`sigma2_imax'

} //end loop 1--over i

local summu_imax=`summu_imax'
local sumsigma2_imax= `sumsigma2_imax'
local dev=(`tsobs'-`summu_imax')/(sqrt(`sumsigma2_imax'))
local pval=1-normal(`dev')


//start calculation of H-L estimate
if "`hlestimate'"=="hlestimate"{
use "`mainfile'",clear

//for consistency in notation with pairsimsens:
local mu=`summu_imax'

//define HL estimate at bound, hat_tau, as (tau_a+tau_b)/2
//tau_a := f(tau_a) <= `mu' & f(tau_a+eps)>`mu' 
//tab_b := f(tau_b) >= `mu' & f(tau_b-eps)<`mu'
//set initial values of tau_l and tau_r
local tau_l=`tau_l0'
local tau_r=`tau_r0'
local j=0
//d is indicator of whether tau_a has been found
local d=0
//first calculate tau_a:
while `d' !=9{
local j=`j'+1

// j indicates run of H-L solver; exit at 1000, indicate failure:
if `j'==1000{
local j=1000 
continue, break
}

if `d'==0 | `d'==-1{
hlarts, intau(`tau_l') resp(`resp') treat(`treat') str(`str')
//hlarts is a subroutine of arsimsens--see below
local ftau_l=`r(hltsobs)'
}

if `d'==0 | `d'==1{
hlarts, intau(`tau_r') resp(`resp') treat(`treat') str(`str')
local ftau_r=`r(hltsobs)'
}

local xw=.5
//(fix weight xw to make slightly faster)
local tau_m=`xw'*`tau_r'+(1-`xw')*`tau_l'
hlarts, intau(`tau_m') resp(`resp') treat(`treat') str(`str')
local ftau_m=`r(hltsobs)'
local tau_mme=`tau_m'-.0000001
hlarts, intau(`tau_mme') resp(`resp') treat(`treat') str(`str')
local ftau_mme=`r(hltsobs)'


if `ftau_m'<= `mu' & `ftau_mme'>`mu'{
local tau_a=`tau_m'
local d=9
}

if `ftau_m'<= `mu' & `ftau_mme'<=`mu'{
local tau_r=`tau_m'
local d=1
}

if `ftau_m'>`mu'{
local tau_l=`tau_m'
local d=-1
}

} //end calculation of tau_a



local tau_l=`tau_l0'
local tau_r=`tau_r0'
local d=0
local j=0
//calculate tau_b:
while `d' !=9{
local j=`j'+1

if `j'==1000{
local j=1000 
continue, break
}

if `d'==0 | `d'==-1{
hlarts, intau(`tau_l') resp(`resp') treat(`treat') str(`str')
local ftau_l=`r(hltsobs)'
} 

if `d'==0 | `d'==1{
hlarts, intau(`tau_r') resp(`resp') treat(`treat') str(`str')
local ftau_r=`r(hltsobs)'
}

local xw=.5
local tau_m=`xw'*`tau_r'+(1-`xw')*`tau_l'
hlarts, intau(`tau_m') resp(`resp') treat(`treat') str(`str')
local ftau_m=`r(hltsobs)'
local tau_mpe=`tau_m'+.0000001
hlarts, intau(`tau_mpe') resp(`resp') treat(`treat') str(`str')
local ftau_mpe=`r(hltsobs)'

if `ftau_m'>= `mu' & `ftau_mpe'<`mu'{
local tau_b=`tau_m'
local d=9
}

if `ftau_m'>= `mu' & `ftau_mpe'>=`mu'{
local tau_l=`tau_m'
local d=1
}

if `ftau_m'<`mu'{
local tau_r=`tau_m'
local d=-1
}


} //end calculation of tau_b

//captured in case tau_a and tau_b not found:
capture local hle=(`tau_a'+`tau_b')/2

//needed b/c round(`x') when " `x' " does not exist returns 0:
if `j' !=1000{
local hlr=round(`hle',.001)
}

} //end hlestimate loop.




use "`mainfile'",clear

ret scalar pval=`pval'
ret scalar dev=`dev'
ret scalar tsobs=`tsobs'
ret scalar gam=`gam'
ret scalar del=`del'
ret scalar summu_imax=`summu_imax'
ret scalar sumvar_imax=`sumsigma2_imax'
if "`hlestimate'"=="hlestimate"{
//captured in case H-L estimate not found:
capture ret scalar hle=`hle'
capture ret scalar hlr=`hlr'
}

} //end quietly

if "`hlestimate'"!="hlestimate"{ 
di  _newline(1) "{result:{ul on}SIMULTANEOUS SENSITIVITY ANALYSIS FOR ALIGNED RANK TEST{ul off}}" _newline(2) ///
"{res:Gamma: `gam'}"  _newline(1) ///
"{res:Delta: `del'}"  _newline(1) ///
"{res:Expectation: `summu_imax'}" _newline(1) ///
"{res:Variance: `sumsigma2_imax'}" _newline(1) ///
"{res:Deviate: `dev'}" _newline(1) ///
"{res:max[p(t>=`tsobs'|u)]: `pval'}"
}

if "`hlestimate'"=="hlestimate"{ 
di  _newline(1) "{result:{ul on}SIMULTANEOUS SENSITIVITY ANALYSIS FOR ALIGNED RANK TEST{ul off}}" _newline(2) ///
"{res:Gamma: `gam'}"  _newline(1) ///
"{res:Delta: `del'}"  _newline(1) ///
"{res:Expectation: `summu_imax'}" _newline(1) ///
"{res:Variance: `sumsigma2_imax'}" _newline(1) ///
"{res:Deviate: `dev'}" _newline(1) ///
"{res:max[p(t>=`tsobs'|u)]: `pval'}" _newline(1)
if `j'!=1000{
di "{res:min[H-L Point Estimate]: `hlr'}"
}
if `j'==1000{
di "{res:Warning: Hodges-Lehmann Estimate Not Found}"
}
}
end


cap program drop elsympolycalc

program define elsympolycalc, rclass
version 9



syntax varname [if] [in]
marksample touse
//local resp="`varlist'"

//find elementary symmetric polynomials 
//e_0(x_1*...x_i*...x_n) ... e_n(x_1*...x_i*...x_n), where x_i is obs i in var `varlist'* 
//via Summation Algorithm
//let e_k=S_k^(i) be denoted below e_kPI

quietly summ `varlist' if `touse'

forvalues N=1/`r(N)'{
 tempname e_0P`N'
 scalar `e_0P`N''=1 //e_0=1
 if `N'>1{
  tempname e_`N'P1
  scalar `e_`N'P1'=0 
 }
}

tempname e_1P1
scalar `e_1P1'=`varlist'[1]
forvalues I=2/`r(N)'{
 forvalues K=1/`r(N)'{
  local Im1=`I'-1
  local Km1=`K'-1
  tempname e_`K'P`I'
  scalar `e_`K'P`I''=`e_`K'P`Im1''+`varlist'[`I']*`e_`Km1'P`Im1''
 } 
}

  

forvalues K=0/`r(N)'{
ret scalar e_`K' = `e_`K'P`r(N)''
}
di "elsympolycalc results returned in r()"

end 





cap program drop hlarts

program define hlarts, rclass
version 9
//hlarts calcuates test stat after proposed tau_a(b) is subtracted from treated responses

syntax [anything], Intau(real) Resp(varname numeric) Treat(varname numeric) Str(varname numeric)

marksample touse 
markout `touse' `treat' `str' `resp'

//for aligned-rank test:
//stratum avg resp:
gen double __hlresp=`resp' if `touse'==1
replace __hlresp=`resp'-`intau' if `treat'==1 & `touse'==1

egen double __hlstavresp=mean(__hlresp) if `touse'==1, by(`str')
//"aligned response":
gen double __hlalresp=__hlresp- __hlstavresp if `touse'==1
//aligned rank:
egen double __hlrespvec=rank(__hlalresp) if `touse'==1
//gkr2000 adjustment:
quietly summ __hlrespvec if `touse'==1
replace __hlrespvec=__hlrespvec/(.5*`r(max)') if `touse'==1

//observed test stat:
egen double __tv=total(__hlrespvec) if `treat'==1 & `touse'==1
egen double __tv2=mode(__tv)
local hltsobs=__tv2
drop __tv __tv2

ret scalar hltsobs=`hltsobs'

drop __hlrespvec __hlalresp __hlresp __hlstavresp

end

